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Abstract 

The periodic standing-wave method for binary inspiral computes the exact numerical solution for 
periodic binary motion with standing gravitational waves, and uses it as an approximation to slow 
binary inspiral with outgoing waves. Important features of this method presented here are: (i) the 
mathematical nature of the "mixed" partial differential equations to be solved, (ii) the meaning 
of standing waves in the method, (iii) computational difficulties, and (iv) the "effective linearity" 
that ultimately justifies the approximation. The method is applied to three dimensional nonlinear 
scalar model problems, and the numerical results are used to demonstrate extraction of the outgoing 
solution from the standing-wave solution, and the role of effective linearity. 

PACS numbers: 04.25.Dm,04.30.Db,04.25.Nx,04.25-g,02.60.Lj 



I. INTRODUCTION 



Background 

The inspiral and merger of a binary pair of compact objects (holes or neutron stars) is one of the most promising 
sources of signals detectable by gravitational wave observatories. For the ground-based detectors LIG0[3|, VIRGOj^, 
GEO6O0Q and TAMAQ, binary merger, especially of intermediate-mass black holes|^ is an exciting possibility; for 
the space-based LISA detector 6, the detection of inspiral/merger of supermassive holes is highly probable, and is 
one of the primary scientific targets. 

The need for theoretical waveforms for the inspiral/merger has driven the effort to find a computational solution 
for the details of the process, but the difficulty of the task has made this problem also a measure of the usefulness of 
numerical computation in general relativity. The hope has been that numerical codes evolving initial data can compute 
the orbital motion using Einstein's equations and, in the case of neutron stars, using hydrodynamical equations. These 
evolution codes would have, as an intrinsic feature, the loss of energy by the binary due to outgoing wave energy, and 
the gradual inspiral due to this loss. 

An important reason for the limited progress on this problem is the matter of timescales. Near a black hole, the 
timescale on which the gravitational field can change is GMjc" , where M is the mass of the hole; for a neutron star the 
timescale is several times longer. The time step in evolution codes is governed by this short timescale. (More precisely, 
the spatial grid near the compact objects must be smaller than GMjt?, and to satisfy the Courant condition, the 
time step must be no larger than 1/c times this grid size.) By contrast, the timescale yjr^ /GM , for orbital motion 
at radius r, is much larger than this, and the timescale for the interesting dynamics, the radiation-reaction driven 
inspiral is much greater yet. The consequence of this incompatibility of timescales is that a very large number of 
time steps is needed in order to see the physics of interest. And computing a large number of time steps is not yet 
possible. Instabilities!^, 01 operating at the short timescale prevent the code from giving useful answers about the 
long timescale. 

The origin of the difficulty suggests its cure: an approximation method that avoids the short timescales. Here we 
describe such a method: a solution for periodic sources and fields. We assume that the compact objects, and their 
fields, rotate with a constant angular velocity (to be denoted Q, below.) This approximation will fail of course, in the 
very latest stages of inspiral merger, when the orbit decays rapidly due to a secular instability or the dynamics of the 
final merger. But that last stage is, by its very nature, rapid; its timescale is only several times that of the shortest 
timescale of the problem. This last stage, then, can plausibly be handled by numerical evolution codes. Indeed, 
evolution codes, especially with perturbation theory handling the final ringdownfl^, are already near to doing this. 
Our goal, then, is a method that can approximate the solution up to the time that numerical evolution codes can take 
over the task that only they can handle. Our approach is not entirely new; it is similar in underlying motivation to 
a method introduced by Detweiler and collaborators [Til IT^ . but our approach is very different in its details and its 
implementation. It is also very closely related to the approach being used by Friedman and his coUaboratorsflSj . 
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Periodic motion and outgoing waves are, of course, impossible in Einstein's theory, both intuitively and mathemat- 
ically^. For this reason, we will solve for standing waves (to be defined and discussed below) in the gravitational 
field. Our periodic standing-wave (PSW) approach, then, will be to find exact (numerical) periodic standing-wave 
solutions of the Einstein field equations and to use these exact solutions as approximations to the physical problem 
of slow inspiral with outgoing waves. 

The most basic ideas behind this periodic standing-wave solution have already been introduced in a previous 
paperfT^, but the implementation there was applied only to two-dimensional models and was limited in other ways; 
in particular, that paper did not discuss the general meaning of standing waves. An general overview of the PSW 
project has also been given[T^. Here we present a more specific discussion, along with numerical results for three- 
dimensional models. This paper is meant to serve as the introduction to the PSW, with subsequent papers presenting 
more detailed information on particular methods, and progress on solving the physical problem. 

Effective linearity and uses of the method 

A key idea in our approach is the relationship of standing waves to outgoing waves. In a linear field theory, a 
definition of standing waves is that they are half the sum of an outgoing solution and ingoing solution. Here, as in 
Ref. we shall call this sum LSIO for linear superposition of (half) ingoing and (half) outgoing solutions. In a 
linear theory, such a superposition is itself a solution. In our nonlinear field model theories, it will turn out that 
— despite strong nonlinearities — this continues to be very nearly true. This effective linearity, the approximate 
equality of the LSIO and a true standing-wave solution, has already been demonstrated for simple two-dimensional 
models, and results for three dimensional nonlinear models will be presented below. More important, the basis for 
effective linearity appears to be robust. This basis lies in the fact that the strong nonlinearities in our model theories 
(and in the physical problem) are confined to the near-field regions around the sources. In these nonlinear near- 
field regions the solution is insensitive to the distant boundary conditions; it is substantially the same for ingoing 
boundary conditions as for outgoing. In this near-field region then, the LSIO will be very nearly a solution despite 
strong nonlinearities, since we are superposing nearly identical solutions. Outside this strong-field near zone the model 
theories, and the physical theory, are nearly linear, so that again the LSIO is a solution. The LSIO will therefore be 
a good approximation to a solution everywhere. 

Below, we shall choose our definitions of standing- wave solution to be close to that of a LSIO, and our approximation 
to a large extent is based on interpreting a standing- wave solution to be approximately a LSIO. In the weak field 
region this LSIO can be deconstructed into outgoing and ingoing pieces and this deconstruction can be continued to 
the strong field source region. (In the source region, the outgoing piece is simply half the solution.) By doubling the 
outgoing piece thus extracted from the standing-wave solution, we thereby arrive at an approximation to the outgoing 
solution for a nonlinear problem. It is in this manner that wc will extract an approximation of an outgoing solution 
from a computed standing-wave solution. 

This is an appropriate place to point out, though not for the last time, the importance of model problems. In 
Einstein's theory there is no obvious meaning to the periodic outgoing solution, so one cannot make statements about 
it, let alone carry out numerical studies. Statements and computations are possible for nonlinear model problems, so 
that tests of effective linearity with such models are crucial. 

The outgoing solutions extracted from our exact periodic solutions can serve two purposes. First, we can use a 
quasistationary sequence of outgoing approximations as a model for the slow physical inspiral. In this approach the 
mass of the system, measured in a weak wave zone far from the orbiting sources, decreases due to the loss of energy 
in outgoing radiation. When we find the system energy as a function of orbit radius, and we compute the outgoing 
radiation, we can infer the rate at which the orbital radius decreases. The difficulty, as with any such quasistationary 
sequence, is how to know that we are comparing the "same" system at different radii. In the case of neutron star's 
the answer is clear; baryon number is an unchanging tag that identifies neutron stars to be the same. For black holes, 
the equivalent tag would be some local mass. The concept of an isolated horizon might give us that local mass. 

The second use for our extracted outgoing solutions is to provide initial data for evolution codes. A spacelike slice 
of our extracted outgoing field will be an excellent approximation to the physical initial data, and should be very 
nearly a solution to the initial value equations. With little change our extracted outgoing initial data, can be made 
into exact (numerical) initial data through the use of York's decomposition 18. 19iJ. 



A periodic binary would have no secular change in its energy, but gravitational waves intuitively rem ove energy. This argument can be 
made mathematically complete using the conservation law for H°'^>^'^, as defined in Misner et al AlS . 
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These two purposes of our solution are not distinct. The natural end point for a quasistationary sequence of PSW 
solutions is the "last orbit," the final stage of motion at almost constant radius. This stage may end due to a secular 
instability, like that of a particle in a black hole spacetime, or due to the imminence of the merger, the formation of 
the final black hole. In either case, this end point must be handled by a numerical evolution code, and in either case, 
the quasistationary sequence will provide ideal initial data for the continuation of the problem by numerical evolution. 

The nature of the mathematical problem 

In the standard approach to computing binary inspiral, initial data are evolved forward in time. In our approach, 
with periodic symmetry imposed, there is no evolution in the usual sense, and there is not the usual concept of initial 
data. Rather we must satisfy boundary conditions at large radius: outgoing, ingoing or standing-wave boundary 
conditions in model problems, and only standing-wave conditions in general relativity. The boundary value problems 
that we must solve differ in two important ways from common boundary value problems. First, our partial differential 
equations (PDEs) are of mixed type. They are of elliptical character in some regions and hyperbolic character in others; 
this will be particularly clear in the model problem to be presented below. We will argue that the mixed character 
causes no fundamental difficulty, and will demonstrate this with the model problems. The mixed character, however, 
does complicate the use of some of the most efficient numerical means of solving bounary value problems. Second, 
we must define what we mean by "standing-wave boundary conditions." Unlike outgoing and ingoing conditions, 
there is no simple local condition corresponding to what we will mean by standing waves in a nonlinear problem. 
We will present two fundamentally different candidates for the standing-wave condition, and here will present results 
of computations with one those of those two choices. (The alternative choice of standing- wave condition is best 
implemented with a special numerical method, and will be presented elsewhere p^.) 

Stepping back from such details, one may be led to ask more fundamental questions about the whole approach. 
Such questions arise especially because the PSW spacetime we compute has some awkward features. Since the exact 
PSW solution contains an infinite amount of gravitational wave radiation, it cannot be expected to be meet the 
asymptotic fiatness conditions of the theorems about the fall-off of fields. But the spacetime is asymptotically fiat in 
that the spacetime curvature decreases with increasing distance from the binary source. Another sign that the PSW 
spacetime has rough edges is that it must not have regular null infinities; Gibbons and Stewart |23| have shown that 
spacetimes with well-behaved Scri-|- and Scri- cannot be periodic. 

It is useful, before diving into details, to clarify what the relationship is between the slightly singular spacetime 
we will be computing, and the physical problem that really interests us. To make this connection we can think of 
the binary system going through several orbits at almost constant radius. A weak wave zone exists at some distance 
from the orbit during this epoch of the motion. The stippled region in Fig. ^ shows the relevant region as part of the 
larger physical spacetime. In this limited region the source motion and the fields are almost periodic, and it is in this 
region only that we hope to approximate the physical fields by the outgoing fields extracted from the computed PSW 
solution. The imperfect asymptotic structure of the PSW spacetime is therefore irrelevant to its physical usefulness. 



time 




FIG. 1: The PSW solution is meant to be an approximation to the physical spacetime only in a limited region. 

In the remainder of this paper we will first present, in Sec^I the mathematical details of a nonlinear model with 
which we clarify many aspects of the PSW approximation. We then discuss, in Sec. lIIII the numerical methods needed 
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to find PSW solutions, especially those aspects of the numerical methods that are idiosyncratic to the special features 
(mixed character, standing-wave boundary conditions, nonlinearities) of our problem. In this section also, results are 
presented of the numerical methods. The results are discussed, and put into the context of the next steps in this 
project [23, in SecHVI 

II. PERIODIC SOLUTIONS, STANDING WAVES, AND MODEL PROBLEMS 

Mixed PDEs and well-posedness 

As stated above, we seek a solution to Einstein's equations in which the sources and the fields rotate rigidly. The 
mathematical statement of this rigid rotation is that there is a helical Killing vector, a Killing vector that is timelike 
close to the sources and spacelike far from the sources. (For more on helical Killing symmetry see Ref. ^2^.) For 
fields in flat spacetime our Killing vector ^ takes the form 

e = dt + nd^ (1) 

in spherical or cylindrical spatial coordinates, and 

i^dt + n{xdy~yd^) (2) 

in Cartesian spatial coordinates. The parameter fi, which must be a constant, can be thought of as the rotation 
rate of the source and fields with respect to an inertial reference frame. For the flat spacetime case, the null surface 
^ • ^ = is a cylinder of radius 1/51 coaxial with the rotation axis. (Here and below we use units in which G ~ c — 1.) 
This cylinder separates the inner region of timelike ^ from the outer region of spacelike ^, as shown in Fig. [21 Since 
this surface, in a sense, represents the points at which the rigidly rotating fields are moving at c, we call this surface 
the "light cylinder," in analogy with pulsar electrodynamics. 

radiative boundary 
]' conditions 




FIG. 2: The "light cylinder" separating the eUiptic and hyperbolic regions of the problem intersects the large spherical surface 
on which numerical boundary conditions are imposed. 

One immediate advantage of the helical symmetry is that it reduces the number of independent variables, thereby 
greatly reducing the computational difhculty of a problem. In our simple flat spacetime models this reduction is most 
easily understood by the fact that helically symmetric scalars cannot depend in an arbitrary way on the spherical 
Minkowski coordinates t, r, 0, 4> but can depend only on t and 4> in the combination ip = (p — Q,t. Thus, in the r, 0, Lp 
system the Killing vector is ^ = dt- These ideas are clarified with a simple flat spacetime model theory for a scalar 
field ^ 

*;a;/3ff"'' + \F = S . (3) 

The term _F'(5',a;") is included to allow for nonlinearity; the constant A adjusts the strength of the nonlinearity. In 
order for a helically symmetric solution 5" to exist, the explicit coordinate dependence of F must be compatible with 
the symmetry. That is, F can have explicit coordinate dependence only on r, 9, and (p. The most natural choice for a 
model would be one in which there is no coordinate dependence, one in which the background spacetime is featureless. 
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We include the possibility of spatial dependence for convenience below. Changing the spatial dependence will help to 
clarify the accuracy of the PSW approximation when nonlinearities are important. 

In the application of the PSW method to holes, an inner boundary condition will be used at a small, approximately 
spherical surface. For simplicity here, however, we use an explicit source term. 



S = S{0 - n/2) [5{^) + - tt)] 



(4) 



representing two points, each of unit scalar charge, in equatorial circular orbits, with radius a and angular velocity f2. 
This source term S obeys the symmetry property that is necessary if a periodic solution is to exist: its Lie derivative 
vanishes along the Killing orbit ^. 

If we are interested only in helically symmetric solutions, then the field equation ||2J) reduces to 



J.2 Qj, 



,5* 

dr 



+ 



1 



d 



'sine 89 



sin 9- 



5* 



1 



sin^ ( 



+ XF{^,r,9,ip)^S{r,9,ip) . 



(5) 



The mixed character of this PDE shows clearly in the coefficient of d^^. The light cylinder is at rsm9 = l/fl 
where this coefficient changes sign. Inside the light cylinder (rsin^ < l/^l) the equation is elliptical; outside it 
is hyperbolic. For outgoing solutions of this equation we impose an outer boundary condition 9,.^' = —dt'^, or 
equivalently dr'^ = fld^p'^, on a spherical surface with a radius large compared to l/fl. As illustrated in Fig.|21 this 
spherical surface is well outside the light cylinder in the equatorial plane, so our boundary conditions are imposed on 
a surface that passes through both the elliptic and hyperbolic regions of the problem. 

Problems with boundary conditions on closed surfaces are common in the case of elli ptic al PDEs. We argue here 
that our boundary value problem with mixed PDEs may be unusual, but is well posed j23l|. Again, a simple model 
problem will help to clarify issues. We set the nonlinearity parameter A to zero in Eq. (O so that we can solve the 
resulting linear equation as an infinite series. If we choose a Dirichlet condition 



4-1 



= 0, 



(6) 



at a finite radius Tmax, then the solution to this linear problem can be written in terms of spherical Bessel functions 
ji,ne as 

* = E E ^'""^'",ilol°^^'"^^'^^ 'K"^»^<) K(m»r>)j,(ml]w) - n,(ml]w)j,(m»r>)] . (7) 



£ m— even 



Here r<(r>) indicates the smaller (greater) of the quantities r, a. Vanishing of the jeijnilr^i^x) denominator means 
that the the "cavity" r < rmax has a resonant mode at frequency fl. In the case that rmax has one of the resonant 
values, the solution to the boundary value problem is not unique. Such values of r are of zero measure, but are dense 
in the set of all r choices. This means that the cavity is always arbitrarily close to a resonance, if sufficiently high 
angular modes are computed. A consequence of this is that a numerical computation does not converge. (Computed 
solutions depend on the computational grid size, and become larger with increasing angular resolution.) The difficulty 
is not just one of computational practice. The boundary value problem is fundamentally ill-posed as a representation 
of fields in an infinite space. There is no meaningful Tmax oo limit of Eq. (|7|). 

Problems with mixed elliptic and hyperbolic regions are of some interest in aerodynamics but there are few 
general results on well-posedness. In those results that do exist, the nature of the boundary conditions plays a pivotal 
role. We have found that this applies to our periodic solutions also. If we replace the Dirichlet conditions of Eq. ^ 
with the Sommerfeld condition 



0, 



(8) 



then the problem is found to be well-posed. This is particularly clear for the linear problem, where the closed form 
solution takes the form 



* = *out + *oxtra ■ 

Here ^'out is the usual "outgoing at infinity" solution 

^out = E E -2imnYi„^iTT/2,0)Ye„^{9,^)ji{mnr<)hf\mnr>) 

£ m— even 



(9) 



(10) 
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and 

*oxtra = Xl -'2imnYf*„^{7r/2,0)Yern{0,ip) jtm jiimnr<)je{mnry) (11) 

£ m— even 

with 

/ h'^^\z) + idh^f^\z)/dz 
^^"^ y jiiz) +idji{z)/dz 

Since the spherical Hankel function has the asymptotic form hg^\z) = (— z)*^^"'"^^e*^[l/z + ©(l/z^)], it follows that 
\jim\ is of order l/rmax- Thus, '5 '^ont, as rmax oo, suggesting that the linear problem is well posed [25l|. 

Numerical results confirm this suggestion. With the boundary condition in Eq. we have encountered no 
fundamental difficulty computing convergent solutions to both the linear and nonlinear versions of Eq. © , and have 
confirmed that solutions do not depend on the particular (large) value chosen for rmax- 




Standing waves: iterative method 



The solutions we will be computing in Einstein's theory, of course, are standing wave solutions, but there are no 
actual "standing wave boundary conditions" analogous to the Sommerfeld condition in Eq. (|SJ) for outgoing waves. It 
is useful, therefore, to explore the meaning of standing-wave solutions with our model nonlinear theories. As pointed 
out in Sec. ^ our paradigm for standing waves is the LSIO of a linear theory, the linear superposition of half ingoing 
and half outgoing solutions. We shall extend this definition of standing wave to nonlinear theories in two ways. The 
first is an extension of the Green function method of Ref. Ts'l, and is called there the TSGF (time symmetric Green 
function) method. For the problem in Eq. |j2Jl this method starts by writing the field equation in the form 



(13) 



Here the operator C['^] depends on but — once 4* is fixed — can be considered to be linearly operating on "if. 
Similarly cjcff depends on but — once 5* is fixed — is to be considered a fixed inhomogeneous term in the equation, 
an effective source term. There is no unique way of putting the field equation into the form of Eq. p3ll for a nonlinear 
model problem, or for general relativity. The quasilinearity of general relativity, and of our nonlinear models, means 
that at least the principal part of C is always unambiguous. There are also some obvious guidelines to follow. In 
particular, C and Coff should become ^'-independent in the weak field limit. 

To iterate for an outgoing solution, for example, one would find an approximate outgoing solution ^Po^t, and then 
would solve 



cm 



](CuV)=^effKJ 



(14) 



for outgoing boundary conditions. The result would be the improved approximation ^'o,i^ to the outgoing waves 



To find standing waves, this method is modified as follows, 
solution. The equation 



An approximation ^'st,j(j is found to the standing-wave 



^[*rtnd](* 



n+l\ 



^offKnd 



(15) 



is then solved with the outgoing boundary conditions of Eq. ((H)) to give ^^out and is next solved with ingoing boundary 
conditions Q — il in Eq. ©] to give ^'"tin^. The new approximation for the standing-wave solution is taken to be 



^n+l ^ 1 ^n+l 
stnd 2 stout 



1 

2 stm 



(16) 



We take the ti — > oo limit of ^'"t^j in Eq. (|16|) to be our computed standing-wave solution. 

We shall call the iterative method just described "direct iteration." This sort of direct iteration is useful in solving 
for the root of an equation x = f{x) only if / is slowly varying. In iteration for ^' the equivalent condition applies 
to C^^ o (Jeff, where is the Green function, the inverse of C for the boundary conditions (ingoing or outgoing) of 
interest. For direct iteration to converge the dependence of o on ^' must be weak and this is the case only 
for nonlinearities of moderate strength. For strong nonlinearities another technique must be used. 

In Newton- Raphson iteration, one uses the iteration 5"" to make linear approximations for £ and Coff. Equation 
(|13|l then takes the form 



/:[^"](^) + (^_^") 



d(JcS 



(17) 
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This equation is linear in ^, and its solution is taken to be the next iteration \I>"+i. By choosing appropriate boundary 
conditions we can use this scheme to iterate an outgoing or ingoing solutions. For our standing-wave solution, we 
follow the same general scheme as with direct iteration. Using the n}^ standing-wave approximation as in Eq. H17(l 
we solve using both in- and outgoing boundary conditions. As in Eq. (|16|l . the n + 1*'' standing- wave approximation 
is taken to be half the sum of the ingoing and outgoing solutions found this way. 

Standing waves: minimization method 

To explain our second, independent way of defining and computing standing-wave solutions, it is best to start with 
the standing-wave solution in the linear model problem. This is simply half the sum of the "outgoing at infinity" 
solution in Eq. and the equivalent ingoing solution. The result is 

^-stnd-JI J2 2mQY;^{7r/2,0)Ye„ie,ip)jiimnr<)neimnr>) . (18) 

£ m— even 

In this solution each multipole has an equal amplitude for ingoing and outgoing amplitudes waves, and one might 
suspect that this property suffices to define standing waves for a nonlinear model. This is not, in fact the case, since 
we could add a multiple of the homogeneous solution ji{mnr)ji(mQa) to the £, m mode without changing the balance 
between ingoing and outgoing. This degree of freedom is equivalent to the degree of freedom inherent in the phase 
between the outgoing and ingoing waves. This extra degree of freedom exists also (though not so transparently) in a 
nonlinear problem. 

To resolve this degree of freedom we can use a generalization of a property that is unique, in the linear case, to 
the correct standing- wave solution Eq. 1)18(1 : In each multipole, the solution is required to have the minimum wave 
amplitude of any solution with balanced ingoing and outgoing waves 26]. 

This method, while very interesting in principle, is difficult in practice to implement in a finite difference boundary 
value approach. One could imagine using a guess for the value of a multipole coefhcient at some outer boundary, 
and then searching for the value that gives the minimum for the amplitudes of the waves in that multipole. In a 
nonlinear problem, the values of each multipole will influence other multipoles, so the search for minimum waves will, 
in principle, be a search in a many dimensional space. 

The real difficulty of this numerical approach is that it uses multipoles as part of the boundary condition. That 
means that multipole coefficients must be extracted. Even in spherical coordinates, the extraction of the multipole 
coefficients involves a weighted sum over all angular grid points. Most important, this sum would not be performed 
as a postprocessing step on a computed solution, but rather would have to be written as a set of equations that would 
form part of the a priori problem to be solved. The set of equations to be solved would then have, in addition to great 
complexity, a boundary-related subset connecting distant grid points. The matrix representation of these equations 
would not have banded structure. In addition to these technical difficulties, the use of spherical coordinates is very 
ill suited to the structure of our source objects, so coordinate patches for the sources would be required. 

For these reasons, we have not attempted to use the minimization criterion in a finite difference code. We have, 
however, implemented this criterion with a spectral approach based on a specially adapted coordinate system. Results 
from this approach are extremely encouraging, but the approach poses new computational challenges, so we are 
continuing to explore two distinct paths: finite difference methods and the iterative definition of standing waves, and a 
spectral/adapted coordinate technique for the minimization criterion. Since the adapted coordinate system necessary 
for the second approach requires a separate development, and is not fundamental to the PSW approximation, we 
confine the present discussion to the first approach, finite difference boundary value problems, with the iterative 
criterion for standing waves. 

III. NUMERICAL IMPLEMENTATION AND RESULTS 

Extraction of an outgoing approximation 

Model problems allow us to test a key idea of the PSW approach, that a good approximation of the outgoing 
solution can be extracted from the computed standing-wave solution 



^'stndcomp = ^ ^ a(rnir)Yim{0,ip) 

eveni m=o,±2,±4.. 



(19) 
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The coefficients Q!(?m(r) are computed from ^'comp, by projection with Y^^. From the reality of 4'stndcomp, the 
coefficients will obey — ctg^ra^ and from the standing- wave symmetry (cosmtp only, no sin m<y9 terms) they will 
also obey aim = a^-m- 

This form of the computed standing-wave solution is compared with a general homogeneous linear (A = 0) standing- 
wave (equal magnitude in- and outgoing waves) solution of Eq. Q, with the symmetry of two equal and opposite 
sources: 



even 



m=0,±2,±4.. 



(20) 



where Ci^m — ^Irn^ from the reality of ^'stndiin- A fitting, in the weak-field zone, of this form of the standing-wave 
multipole to the computed function ai„i{r) gives the value of C'lm- 

By viewing the linear solution as half-ingoing and half-outgoing we define the extracted outgoing solution to be 



^'exout = ^ ^ Yirn{9,ip)Clrnhf\mU.r) 

even£ m=o,±2,±4.. 



(21) 



Since this extracted solution was fitted to the computed solution assuming only that linearity applied, it will be a 
good approximation except in the strong-field region. In the problems of interest, the strong-fields should be confined 
to a region near the sources. In those regions, small compared to a wavelength, the field will essentially be that of 
a static source, and will be insensitive to the distant radiative boundary conditions. As pointed out in Sec. the 
solutions in this region will be essentially the same for the ingoing, outgoing, and standing-wave problem. In this 
inner region then, we take our extracted outgoing solution simply to be the computed standing-wave solution, so that 



X] YirnCirnh 
^stndcomp 



(1) 



weak field outer region 
strong field inner region 



(22) 



The boundary between a strong field inner region and weak field outer region would ideally be a closed surface 
surrounding each of the source regions. This is easily implemented with the adapted coordinates to be introduced 
in a subsequent paper. Here, for simplicity, we take the boundary to be a spherical surface around the origin. In 
order for the extracted solution to be smooth at this boundary, we use a blending of the inner and outer solution in 
a transition region extending between radii riow and fhigh and, in this region, we take 



Here 



/3(r) 



(3{r) = 3 



[1 - /3(r)] ^stndcomp ■ 



r - rio 



n 2 



^high ^low 



- 2 



r - Ho 



^high ^low 



(23) 



(24) 



so that f3{r) goes from at r = riow to unity at r = ^high and has a vanishing r-derivative at both ends. 

In the case of our typical choice SI = 0.3, the value of rjow is chosen to be r = 1.3 a the value at which the static 
and standing- wave solutions of the linear problem differ by 10%. This value should decrease with increasing fi, but 
it must be larger than the orbital radius r = a, so we choose it to be 



How = a[l + 0.3(0.3/17)] . 
In order to have a moderately thin transition region we somewhat arbitrarily take 

r-high = fl[l + 0.6(0.3/17)] . 



(25) 



(26) 



For the numerical results reported below, the extraction details of Eqs. (|19|I - H26|I are used, and extraction is carried 
out using the i — 0,2,4 multipoles. 



Choice of model 



To verify and demonstrate several innovative features (well-posed mixed boundary value problem, standing waves, 
effective linearity) of the PSW approximation, we use the nonlinear scalar model of Eq. jSjl, with the 6 function 
sources given by Eq. Q . We make the simplifying assumption that the nonlinear function F in Eq. depends only 
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on ^, not on its derivatives. 
We replace Eq. by 



From this an obvious simplification follows for the iteration method of Eqs. H13|l - H16|l . 



with C taken to be 



j,2 Qj. K Qj. 



1 d f . d_ 

sine 80 961 



1 



sin^ I 



The effective source term includes both the true point source and the nonlinear term 

5{r 



-S{e - 7r/2) [Siip) + Siip - tt)] - XF . 



(27) 



(28) 



(29) 



Our choice of the nonlinearity function F is 



F 



^5 



\I/4 



(30) 



(We will comment below on the difference between this choice and that made in previous work, including previous 
versions of this paper.) Here \l/o is a second nonlinearity parameter (A being the first). We shall choose "ifo to be less 
than unity; in the numerical results to be presented, ^'o is taken to be 0.15. 

To understand the effect of this nonlinearity, let R denote the distance from one of the point sources. Very near a 
source point, at very small R, where the field is strong, F has the limit F so that the solution of Eqs. H27 |l - H30(l 

approximately has the Yukawa form 



-A R 



AnR 



near source pt . 



(31) 



At some distance from the source — call it i?iin — the field ^ becomes smaller than ^^o, and F can be approximated 
as Vl/^. Since \1/ itself is less than '^q, and hence less than unity, this nonlinearity is small enough to be considered 
a perturbative correction. 

If the transition at i?iin takes place well inside the near zone of the problem, then the effect of the nonlinearity 
can be understood as follows: Near a source point the solution has the form of a unit strength Yukawa potential. At 
distance i?iin the effect of the XF term is turned off and the solution becomes a simple Coulomb potential. The source 
strength for this Coulomb field, though, will be less than unity. Due to the exp {—^/—XR) Yukawa factor, the source 
strength decreases in the region from _R = to i? = i?iin, and the effect of the nonlinearity is to reduce the effective 
source strength by a factor of order exp (— V— Ai?iin). Since this transition takes place well within the near zone, 
it should be this reduced source strength that is responsible for generating radiation. The effect of the nonlinearity 
on radiation, then, will be the same reduction factor exp {—^/—XRun), and we can easily estimate the size of this 
nonlinear effect. One estimate can be found by solving 



inRii, 



= *r 



for i?iin, and using this value of i?iin in the expression exp (— V^-^^iin) for the reduction factor, 
follows by solving the spherically symmetric static nonlinear problem for a unit strength source 



(32) 



Another estimate 



]_d_ 

J.2 Qj, 



IT 1 +^ 

or 



(33) 



(Here the right hand side is the unit 5 function at the origin.) For this solution the ratio is found of the large-r 
monopole moment to the small-r monopole moment, and this ratio is taken as the reduction factor. Since these 
methods for the reduction factor ignore the nonlinear interaction between the two source points, and since they 
assume that all the wave generation occurs far outside i?iin, they can only be considered an approximation for the 
nonlinear reduction effect on the wave amplitude. We shall see, however, that these estimates are accurate enough to 
be taken as a good heuristic explanation of the role of the nonlinearity. 

In previous work, a form of the nonlinearity was used that was different from that in Eq. I|3()() . To give that previous 
form we first defined the distance (i?-) from the source point on the x axis at a; = a (x = —a) to be given by 



i?± = (r sin 9 cos (p =F a)^ 



sin^t 



COS^ (fi - 



cos^6 



(34) 
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We then introduced the distance variable 

X=/R7i?~. (35) 

At either of the source points x ^ 0, and far from the sources x ^ Ii^ terms of x, the form of the nonhnearity 
previously used|22| is 

^P-e.= (£)"e("-/^)^. (36) 

The x-dependent prefactor (x/na)"e*^"~'^/°' was included so that we could force the nonhnearity to be concentrated 
near x = na. By choosing n to be 5 or 10 we could, in this way, have strong nonhnearity in the wave zone, and 
we could numerically demonstrate the failure of effective nonhnearity. The x-dependent prefactor, however, makes it 
difRcult to find a numerical solution that is physically meaningful. 

The prefactor is a difhculty because the solution near the source can have either the Yukawa form exp {—^J—\R), 
or the "anti- Yukawa" form exp (+\/— Ai?). If there is any of the latter included in the solution, then the field gets 
larger at larger distances from the sources, so the strong nonhnearity is never suppressed, the \F term continues 
to approximate and the sum of the Yukawa and anti- Yukawa forms continues to be a valid solution. But if 
the anti- Yukawa part is present, the solution cannot meet the fall-off conditions at an outer boundary at large r. 
Without the prefactor, then, the outer boundary conditions act to suppress the anti- Yukawa part of a solution. With 
the prefactor present, however, the nonhnearity can be turned off by the fall off of exp(— x/a), even if the solution 
contains an anti- Yukawa part close to the sources. The prefactor, in effect, shields the inner region from the influence 
of the outer boundary conditions. When the prefactor is included in the nonhnearity, the solution in the inner region 
will be a somewhat unpredictable mixture of Yukawa and anti- Yukawa parts that is sensitive to grid spacing. 

The choice made for the dependence in Eq. H30|) . rather than that in Eq. I|3t)|) . is motivated by the fact that 
p falls off rather slowly in the weak wave zone. Changing the form of F to ^1^^/(1 -I- '^^) cures this slow fall 

off, but imposes a very sharp cutoff near the sources, one that is too sharp for our relatively coarse computational 
grid. By taking F proportional to 5*^/ (^Pq -|- ^f''), with a fairly small value of 'So, the falloff of F is smoothed out and 
moved to a larger distances from the source. 



Numerical methods 

Since C is independent of ^ we can (as in Ref. ^^) compute once and for all the inverses of C, i.e. , the Green 
functions corresponding to specific boundary conditions. In this way, we can compute and , the Green 
functions for outgoing and for ingoing boundary conditions. The direct iterative method of Eqs. (|14|l - l|15|l . then 
amounts to 

*ouV = Cut ('^cff [*^ut]) = (^cff[*r„]) (37) 



(38) 



Since C has no 5" dependence, the basic Newton-Raphson iteration simplifies to 



C- 



= creff[*"] - 



da, 



cff 



(39) 



This Newton-Raphson approach can be applied to find outgoing, ingoing and standing-wave solutions analogous to 
Eqs. and 

Each iteration of Eqs. (|37|1 . (|38|l or (|39() is equivalent to the solution of a large set of linear equations. Such systems 
are most typically encountered for elliptic boundary value problems, and are typically solved most efficiently with 
relaxation methods, or related methods (e.g. , multigrid) based on the geometry of the problem. Such methods start 
with an approximate set of values for each of the unknowns at every point of the numerical grid. At each point 
the solution is then recalculated on the basis of the values at nearby grid points. This method sweeps through all 
the points of the grid and is iterated until an error criterion is met. Such a method must be compatible with the 
domain of dependence for the points of the grid. For an elliptic PDE, for example, the values of unknowns are 
updated at a central point of a set of grid points. For a hyperbolic PDE, on the other hand, the field computation, 
or updating, must be done only at a point in the "future" of those grid points being used. For a mixed boundary 
value problem a relaxation method has special difficulties, especially at the interface between elliptic and hyperbolic 
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regions. Nevertheless, relaxation methods have been successfully applied to mixed PDEs in transonic aerodynamics, 
first by Murman and Cole[2^. The slow convergence of this method at the interface (the "sonic surface" in transonic 
aerodynamics) can be improved with special techniques that may need to be specific to the Droblem|29l|. 

We are presently investigating relaxation and other numerical methods (e.g. , decomposing the grid into regions and 
applying different techniques, preconditioners, etc.) for large grids and many variables. For our three-dimensional 
scalar problem illustrated here, however, we have been able to use a more-or-less straightforward method of inverting 
the matrix for the finite difference equations. 

In one approach to finding an iterative solution, we use matrix inversion at each step of the direct iteration of 
Eqs. H37() . (|38|l . and we take advantage of the fact that is rotationally symmetric (i.e., it is translationally 
symmetric with respect to (^), and we work with the Fourier components 4'J^(r, 0)e*'"'^ of the iterative solution. At 
each step of iteration we project out the Fourier components of the effective source. Due to the nonlinearity in the 
effective source, the Fourier modes of ^'^(r, 6')e""'^ mix in this step, but is rotationally symmetric so the Fourier 
modes do not mix in the step of solving for This method takes advantage of the efficiency of a fast Fourier 

transform (FFT) and reduces the RAM needed to little more than that for a two-dimensional r, d grid. This method, 
therefore, allows a rather fine grid in r and 0. 

We have used this efficient FFT method extensively, but direct iteration has the drawback already cited following 
Ea. (|16|l : it is limited in the strength of the nonlinearity it can handle. Direct iteration will not converge for very strong 
nonlinearity. The iterative Newton- Raphson method of Eq. 139f) . on the other hand, docs almost always converge once 
one has a solution sufficiently close to the correct solution. The operator on the left in Eq. H39|) . however, contains the 
previous iteration which is not symmetric in ip, so that the FFT method cannot be used with Newton- Raphson 
iteration. This has meant that a relatively coarse grid, or large RAM, had to be used. We have not yet implemented 
a parallelizable method for solving the iteration steps, and have been restricting most runs to 8GB. 

It is worth mentioning an interesting hybrid method that we have explored. The problem in Eqs. H27|) - (|30|l . outside 
the point sources, can be written as 

(^ + ^)* = ^*o^4^- (40) 

The nonlinearity on the right is never large; it is small both for > 'So near the sources, and for < far from the 
sources. The weakness of the formal nonlinearity suggests that a solution may converge with direct iteration even for 
large nonlinear effects. The operator (£ -|- A), furthermore, is rotationally symmetric, so the FFT method can be used. 
The method, however, turns out to have a serious flaw. Where ^' is small, the left-hand side of Eq. H40|) is dominated 
by A^*, which is very nearly equal to the right-hand side. In the analogy we gave, following Ea. (|16|l . to the iterative 
solution for a root of a: = /(x), this is equivalent to / having a derivative very close to unity at the root. It would 
appear that this difficulty could be avoided by iterating Eq. I|4(J|) near the source, where the nonlinearity is strong, 
and iterating the standard form of the problem in the weak field region. Numerical experiments with this approach 
have been inconclusive. Since we do not intend to use a single-patch spherical coordinate sysem in the future, we 
have not examined this hybrid method exhaustively. 

In practice we have used the following eclectic approach to find solutions: (i) For linear models, for which no 
iteration is required, we have taken advantage of the RAM-reduction of the FFT method, (ii) For strongly nonlinear 
models we have used Newton-Raphson iteration on a three-dimensional (non-FFT) grid, and have used continuation 
(i.e. , "ramping up") both in A and in fi. Despite RAM limitations, we have been able to confirm that the solutions are 
second-order convergent, (iii) For —A less than around 2, it has been been possible to find solutions with the direct- 
iteration, FFT method. These solutions have been compared with the corresponding solutions from the non-FFT, 
Newton-Raphson method, and have been found to agree within the numerical uncertainty in the solutions. 

In applying the iteration methods, and looking for convergence, we have used two error measures. One, eitor = 
xj/" _ is the difference at a grid point between the computed value at a grid point, and the value computed 

at the previous iteration. The second error measure, Csoin is the value of (£(^'") — creft[5'"]) at a grid point. In our 
FFT computations the criterion for convergence was to have the rms value of Citcr (averaged over the entire grid) fall 
below 1 X 10~^. The value of esoin was also monitored in the FFT computations and was found not to be larger than 
1 X 10^^ at any grid point, and to have an rms value typically around 1 x 10~^. A much more stringent requirement 
for convergence was used in the Newton-Raphson computations: the rms value of esoin bad to fall below 5 x 10^^^ for 
the solution to be acceptable. 

Numerical results 

We first illustrate the fundamental concept of the PSW method with various solutions of Eq. ©, with the nonlin- 
earity given in Eq. 13U|) . Figure |21 shows solutions in the equatorial (Q = 7r/2) plane; the amplitude of the field ^' is 
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FIG. 3: The ^ field for two rotating point sources in the equatorial plane. The fields shown are nonlinear solutions of Eqs. 
and II3UII . with aQ = 0.3 and A = —1. For clarity, the (^-average is removed at each radius. Parts (a) and (b) of the figure 
show, respectively, the nonlinear outgoing and ingoing solutions. Part (c) is the standing-wave solution, and part (d) is the 
outgoing solution extracted from it. The vertical scale gives field strength (arbitrary units) and the horizontal coordinates are 
corotating Cartesian coordinates in units of a, the distance of a source from the rotation axis. 

plotted as a function of corotating Cartesian coordinates x — rcosip, and y — rsiriLp. The source points are on the x 
axis at a; = ±a, and the outer boundary is at r = 30a. For all four plots, aft — 0.3, A = — 1 and ^'o — 0.15. The results 
plotted are those from direct iteration with the FFT method, for a computational grid using 361 radial divisions, 16 
divisions in 9, and 32 Fourier modes. For all models, the computed results are dominated by the monopole, so for 
clarity in the figures the (/3-average of the solution has been subtracted at every radius. It is worth noting that this 
procedure not only removes the monopole (the £ = part of the solution), but also removes the m = part of the 
quadrupole, etc. 

The plot in part (a) of Fig. O shows the outgoing solution (solution for outgoing boundary conditions); the plot 
in part (b) shows the corresponding ingoing solution. The plot in part (c) is the computed standing-wave solution 
for the same problem parameters. (Note: this nonlinear standing-wave solution is not half the superposition of the 
outgoing and ingoing solutions. Rather, it is the nonlinear field equation solved with the standing-wave definition 
discussed in Sec. E]) Part (d) shows the key idea of the PSW approximation, the outgoing solution extracted from 
the standing- wave solution, by the extraction method described in Eqs. (|19|) - (|26|l . When the PSW method is used in 
general relativity, it will be possible only to compute the standing-wave solution; the extracted outgoing solution will 
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TABLE I: The reduction factor for outgoing waves due to the nonlinearity. For all cases, 'i'o — 0.15. The second column refers 
to Eq. I32II . This equation is solved for Run- Estimate 1 uses this value of i?iin in exp (— V— Ai?iin). Estimate 2 is the reduction 
factor found from a numerical solution of Eq. 13311 . The last column gives the results from Newton-Raphson computation with 
aQ = 0.3, with the outer boundary at r = 30a, and with a. r,9,(p grid of 120 x 20 x 32. 



A 


Estimate 1 


Estimate 2 


Computation 


-1 


69% 


87% 


78% 


-2 


62% 


73% 


68% 


-5 


53% 


65% 


55% 


-10 


46% 


54% 


47% 


-25 


37% 


41% 


35% 



TABLE 11: Convergence of finite difference computations. Nonlinear outgoing solutions are computed with five different grid 
resolutions for A = — 10, ^'o = 0.15, afi = 0.3 and outer boundary at r = 30a. An L2 norm is computed for the difference 
between the solution for grid k and grid k + 1. This is reported as the "Error" for grid k. 



k 




Error 


1 


90x10x16 


2.71 E-5 


2 


120x14x22 


1.60 E-5 


3 


150x16x26 


8.68 E-6 


4 


180x20x32 


5.22 E-6 


5 


210x24x38 





represent the approximation to the physical, outgoing solution. 

Table gives quantitative results for strongly nonlinear outgoing waves. In that table, values are given for the 
reduction factor due to the nonlinearity. As explained following Eq. H31|) , this is the factor by which the nonlinearity 
decreases the amplitude of the waves. (For the same aU, and source strength, the amplitude of outgoing waves for the 
linear problem A = problem is compared to the amplitude for a problem with X ^ Q.) The fact that the reduction 
factors are significantly different from 100% shows that we are able to compute models in which nonlinear effects are 
strong. In the table the computed reduction is compared with estimates from heuristic models of Eqs. (|31|l - H33f) in 
which \E' is taken to have a Yukawa form very near the source, and a Coulomb form further out, but well within the 
near zone. The agreement of the computations with the estimates is strong evidence that the heuristic model captures 
much of the nature of the nonlinear effect. 

Table ^1 gives information on the numerical errors of the most computationally intensive solution type: that for 
strongly nonlinear waves computed via Newton-Raphson iteration. A single physical model (A = —10, ^'o = 0.15, 
ail = 0.3, outer boundary at r = 30a) is computed on five different grids. As a measure of the truncation error for 
grid fc, the L2 difference (the square root of the average square difference) is found between the results for grid k and 
for grid fc + 1 . This is listed in Table ITU as the error in grid k. These results, especially for the finest three grids, 
suggest quadratic convergence of the numerical process. 

The crux of the PSW method is that a good approximation to a nonlinear outgoing solution can be extracted from 
a standing wave nonlinear solution. Examples of this are given in the next two figures, the central numerical results 
in this paper. 

Figure^ shows results for computations of linear (A = 0) and nonlinear (A = —2 and -10) models for 5" along 
the 9 = 7r/2, ip ~ Q lines. All models used rotation rate afi = 0.3 and nonlinearity parameter = 0.15 and the 
180 x 20 x 32 grid with an outer boundary dX r — 30a. The nonlinear models were solved through Newton-Raphson 
iteration with continuation in A. The figure shows, as continuous lines, the computed solutions for outgoing waves 
with A = 0,-2, —10. Included in the figure also are the A = —2, —10 results for the approximate outgoing waves 
extracted from the nonlinear standing wave solutions by the method of Eq. H23|l . The difference between these 
outgoing approximations and the true outgoing waves is so small that the approximation results are given as discrete 
data points to aid in visualization. 

Figure |S1 shows the small-radius portions of the same models as those in Fig. 21 (The A = curve is nearly 
indistinguishable from that for A = —2, and is omitted from the figure.) The radial 6 — 7r/2, = Q line along which 
the results are presented, goes through the source point at r = a, so Fig. [31 shows the computed solution in the 
neighborhood of the source. The Yukawa-like effect of the nonlinearity near the source is evident in more rapid fall-off 
of the A = — 10 model away from the source point. 

The results in Figs. [3 and 01 are graphical evidence for the accuracy of the PSW method; the outgoing waveforms 
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FIG. 4: Extracted outgoing nonlinear waves vs. true outgoing nonlinear waves. For A = 0, —2, —10, \l/o ~ 0.15, af2 = 0.3, with 
a 180 X 20 X 32 grid. The field ^ is shown as a function of r along a radial line through the source point, i.e. , along the 9 = 7r/2, 
if = line. Continuous curves show computational results for outgoing waves. Discrete points, for the nonlinear models, show 
the approximate outgoing waves extracted from standing wave solutions. 




-2 
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FIG. 5: The same models as in Fig.|l] but in the region of the sources. As in Fig. |1] continuous curves show the computations 
of the true nonlinear outgoing waves, and discrete points show the outgoing wave approximation extracted from the nonlinear 
standing wave solution. 

extracted from the nonlinear standing wave solution are excellent approximations to the true outgoing waves both 
near the sources and in the wave zone. The agreement in the intermediate zone (not shown in the figures) is equally 
impressive. A quantitative measure of the agreement is the "L2" difference of the outgoing wave and the extracted 
outgoing approximation. This measure is the square root of the average (over all grid points) of the squared difference 
between the true and the extracted outgoing solutions. For A = —10 this L2 difference is 8.7 x 10^^ and is of the 
same order as the error in Table Hll for the 180 x 20 x 32 grid being used. Since the numerical uncertainties are of the 
same order as the difference between the true and the extracted outgoing waves we cannot claim to have computed 
any meaningful inaccuracy in the PSW approximation. 

This is unfortunate. In presenting numerical results it would be useful to demonstrate that the PSW approximation 
is, after all, an approximation by showing a model in which the extracted outgoing solution is significantly different 
from the true outgoing solution. Our inability to do this is related to limitations on numerical solutions. Our arguments 
for effective linearity show that the PSW approximation should fail only if the region of significant nonlinearity overlaps 
the wave zone. For this reason we used the x-dependent prefactor of Eq. H3t)|l in an earlier version of the present paper 
to allow us to force the nonlinearity to be concentrated in the wave zone. Although that technique did allow us to 
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induce significant errors in the PSW approximation, we have explained, following Eq. H36() . why the solutions for the 
models with the prefactor have undesirable features. If no unnatural x-dependence is explicitly injected in the source, 
the way in which effective nonlinearity can be made to fail is for a^l to approach unity, i.e. , for the source points to 
move very relativistically (a case in which the PSW would be expected to fail for binary inspiral). Unfortunately, we 
have not been able to find convergent solutions for large aQ. Presumably this is due to the fact that large ail means 
fields with sharp gradients, too sharp to be handled by our necessarily coarse grids. 

IV. CONCLUSIONS 

We have given here the foundations of the PSW method based on the extraction of an outgoing solution from a 
computed standing-wave solution. We have also given the details of the extraction calculation. 

The results provided for convergent nonlinear models are "proof" by example that there is no fundamental mathe- 
matical problem of well-posedness of the mixed PDE problem, with radiative boundary conditions on a sphere that 
is in both the elliptic and hyperbolic regions of the problem. We have, furthermore, presented limited numerical 
evidence for the validity of the PSW method, i.e. , that the extracted outgoing solution is a good approximation to 
the true nonlinear outgoing solution. This evidence helps make the case for the application of the method to the 
general relativistic problem, in which only the standing-wave solution will be computable, and the extracted solution 
will be taken as the approximation to the physical problem. 

The numerical studies have also taught a lesson about the limitations of the relatively straightforward numerical 
method used here, matrix inversion of the finite difference equations in spherical coordinates. We have found that 
this method is limited by the coarse grid that can be used for the finite-differencing. We could, in principle, use a 
software engineering approach to increase the range of nonlinearity and rotation rate that can be handled. But the 
methods used here, spherical coordinates and delta function sources, are meant only to provide a relatively simple 
context for establishing the foundations for more advanced approaches. 

In a paper now in preparation j23, we will present an important step forward in dealing with PSW problems, a 
coordinate system that conforms to the geometry near the sources and far from the source asymptotically goes to 
spherical polar coordinates, the coordinates best suited to the description of the waves. One advantage of this method 
is that it allows us very simply to put in details of the sources as inner boundary conditions rather than point sources. 
In addition, the new coordinates turn out to be very well suited to a spectral method that has shown remarkable 
computational efficiency, but that poses new computational problems. Computations using an adapted coordinate 
system have already been carried out for the three-dimensional nonlinear scalar problem with both the finite difference 
and spectral formulation, and for linearized general relativity using the finite difference formulation. Since the details 
of adapted coordinates, especially with the unusual spectral method, are not directly related to the foundations of 
the PSW method, those details are appropriate to a separate paper. 

A very different approach to better numerics is to use relaxation methods, already mentioned in Scc. lIIII In view of 
the large number of uncertainties about their application, we have started on a basic study of relaxation methods in 
mixed PDE systems in PSW-type problems, but will continue to explore a number of different numerical approaches 
to the PSW problem. 
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